# Trade Agreement Politicization Analysis
# Author: Scott Hamilton
# Date: October 2025
# Description: Analysis of politicization patterns across EU trade agreements

# Load required libraries
library(dplyr)
library(tidyverse)
library(ggplot2)
library(patchwork)
library(ggcorrplot)
library(broom)
library(olsrr)
library(extrafont)
library(showtext)
library(car)
library(stargazer)
library(nnet)
library(kableExtra)
library(glmnet)


# Set up fonts for publication-quality figures
font_add("cormorant", "fonts/Cormorant-Regular.ttf")
showtext_auto()
font_add_google("Cormorant", "Cormorant")
theme_set(theme_minimal(base_family = "Cormorant", base_size = 14))

# =============================================================================
# DATA PREPARATION
# =============================================================================

# Load and prepare main datasets
salience <- read.csv2(file="data/search_interest_wide.csv")
scrutiny <- read.csv2(file="data/parliamentary_scrutiny_wide.csv") 
media <- read.csv2(file="data/media_salience_wide.csv")

# Standardize column names
colnames(salience) <- c("member_state", "ttip", "ceta", "mercosur", "ukraine", "vietnam", "south_korea", "singapore", "uk", "mexico")
colnames(scrutiny) <- c("member_state", "ttip", "ceta", "mercosur", "ukraine", "vietnam", "south_korea", "singapore", "uk", "mexico")
colnames(media) <- c("member_state", "ttip", "ceta", "mercosur", "ukraine", "vietnam", "south_korea", "singapore", "uk", "mexico")

# Convert to long format and combine
salience_long <- gather(salience, agreement, salience, ttip:mexico, factor_key=FALSE)
scrutiny_long <- gather(scrutiny, agreement, scrutiny, ttip:mexico, factor_key=FALSE)
media_long <- gather(media, agreement, media, ttip:mexico, factor_key=FALSE)

only_salience <- subset(salience_long, select=c("salience"))
only_scrutiny <- subset(scrutiny_long, select=c("scrutiny"))
df <- cbind(media_long, only_scrutiny, only_salience)
df <- df %>% rename(Agreement = agreement, search = salience)

# Create composite measures
df$dyad <- paste(df$member_state, df$Agreement)
df <- transform(df, salience = ((search * 2) + media)/3)
df <- transform(df, contestation = ((scrutiny * 2) + media)/3)
df <- transform(df, politicization = pmin(salience, contestation))

# Define custom labels and shapes for consistency
custom_labels <- c(
  "ceta" = "CETA", "mercosur" = "EU-Mercosur", "mexico" = "EU-Mexico",
  "singapore" = "EU-Singapore", "south_korea" = "KorEU", "ttip" = "TTIP",
  "uk" = "UK Withdrawal", "ukraine" = "EU-Ukraine", "vietnam" = "EU-Vietnam"
)

agreement_shapes <- c("ttip" = 1, "ceta" = 2, "mercosur" = 3, "ukraine" = 4, 
                      "vietnam" = 20, "south_korea" = 23, "singapore" = 15, 
                      "uk" = 18, "mexico" = 8)

# =============================================================================
# CHAPTER 1: DESCRIPTIVE ANALYSIS (Figures 1.1-1.9)
# =============================================================================

# Figure 1.1: Distribution of Politicization
###!!!FIGURE 1.1!!!###
ggplot(df, aes(x = politicization, fill = Agreement, linetype = Agreement)) + 
  geom_histogram(binwidth = 10, alpha = 0.8, position = "dodge", colour = "black") +
  scale_fill_grey(start = 0.1, end = 0.9, labels = custom_labels) +
  scale_linetype_manual(
    values = c("solid", "dashed", "dotted", "dotdash", "longdash", 
               "twodash", "solid", "dashed", "dotted"), 
    labels = custom_labels
  ) +
  theme_bw(base_family = "Cormorant", base_size = 14) +
  labs(x = "Politicization", y = "Count", fill = NULL, linetype = NULL) +
  scale_x_continuous(breaks = seq(0, 100, by = 10)) +
  theme(
    legend.position = "top", legend.title = element_blank(),
    legend.justification = "center", legend.box.just = "center",
    legend.margin = margin(6, 6, 6, 6),
    panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
    panel.background = element_blank()
  ) +
  guides(fill = guide_legend(ncol = 4), linetype = guide_legend(ncol = 5))

# Figure 1.2: Politicization by Agreement
###!!!FIGURE 1.2!!!###
ggplot(df, aes(x = Agreement, y = politicization)) + 
  geom_boxplot() +
  theme_minimal(base_family = "Cormorant", base_size = 14) +
  theme(axis.text.x = element_text(angle = 25, vjust = 0.5, hjust = 1)) +
  labs(x = "", y = "Politicization") +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank())

# Figure 1.3: Politicization by Member State
###!!!FIGURE 1.3!!!###
ggplot(df, aes(x = member_state, y = politicization)) + 
  geom_boxplot() + 
  theme_minimal(base_family = "Cormorant", base_size = 14) +
  labs(x = "", y = "Politicization") +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank())

# Figure 1.4: Correlation Matrix
df_cor_pol <- subset(df, select=c("search", "media", "scrutiny", "salience", "contestation", "politicization"))
df_cor_pol_plot <- round(cor(df_cor_pol, method="pearson"), 6)
colnames(df_cor_pol_plot) <- tools::toTitleCase(colnames(df_cor_pol_plot))
rownames(df_cor_pol_plot) <- tools::toTitleCase(rownames(df_cor_pol_plot))

###!!!FIGURE 1.4!!!###
ggcorrplot(df_cor_pol_plot, 
           hc.order = FALSE, type = "upper", outline.col = "black",
           lab = TRUE, lab_size = 3, lab_col = "black",
           colors = c("black", "gray", "white")) +
  theme_minimal(base_family = "Cormorant", base_size = 14) +
  labs(x = NULL, y = NULL) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), legend.position = "right")

# =============================================================================
# TRAJECTORY ANALYSIS (Figures 1.5-1.9)
# =============================================================================

# Create trajectory datasets
data <- df
data$ID <- seq.int(nrow(data))
names(data)[5] <- "search"

`%nin%` = Negate(`%in%`)

# Define trajectories
unclassified <- subset(data, search == media & media == scrutiny)
unclassified$trajectory <- "unclassified"

df_csd <- subset(data, search >= media & media >= scrutiny)
df_csd <- subset(df_csd, df_csd$ID %nin% unclassified$ID)
df_csd$trajectory <- "Societal Demand"

df_es <- subset(data, search <= media & media <= scrutiny)
df_es <- subset(df_es, df_es$ID %nin% unclassified$ID)
df_es$trajectory <- "Elite Supply"

df_pp <- subset(data, search >= media & media <= scrutiny)
df_pp <- subset(df_pp, df_pp$ID %nin% unclassified$ID)
df_pp <- subset(df_pp, df_pp$ID %nin% df_es$ID)
df_pp <- subset(df_pp, df_pp$ID %nin% df_csd$ID)
df_pp$trajectory <- "Partisan Representation"

df_mz <- subset(data, search <= media & media >= scrutiny)
df_mz <- subset(df_mz, df_mz$ID %nin% unclassified$ID)
df_mz <- subset(df_mz, df_mz$ID %nin% df_es$ID)
df_mz <- subset(df_mz, df_mz$ID %nin% df_csd$ID)
df_mz$trajectory <- "Mediatization"

# Create trajectory plotting function
create_trajectory_figure <- function(data, title) {
  data_long <- gather(data, indicator, value, media:search, factor_key=TRUE)
  data_long$indicator <- factor(data_long$indicator, 
                               levels = c("search", "media", "scrutiny"),
                               labels = c("Search", "Media", "Scrutiny"))
  
  ggplot(data = data_long, aes(x = indicator, y = value, group = dyad)) +
    ggtitle(title) +
    geom_line(aes(linetype = Agreement)) +
    geom_point(aes(shape = Agreement)) +
    scale_shape_manual(values = c(1:28)) +
    scale_linetype_manual(values = c("solid", "dashed", "dotted", "dotdash", 
                                     "longdash", "twodash", "solid", "dashed", "dotted")) +
    theme_minimal(base_family = "Cormorant", base_size = 14) +
    labs(x = "Indicator", y = "Value", linetype = "Agreement", shape = "Agreement") +
    theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
          legend.position = "right")
}

# Figure 1.5: Societal Demand
###!!!FIGURE 1.5!!!###
create_trajectory_figure(df_csd, "Societal Demand")

# Figure 1.6: Elite Supply  
###!!!FIGURE 1.6!!!###
create_trajectory_figure(df_es, "Elite Supply")

# Figure 1.7: Partisan Representation
###!!!FIGURE 1.7!!!###
create_trajectory_figure(df_pp, "Partisan Representation")

# Figure 1.8: Mediatization
df_mz_long_copy <- gather(df_mz, indicator, value, media:search, factor_key=TRUE)
df_mz_long_copy$indicator <- factor(df_mz_long_copy$indicator, 
                                   levels = c("search", "media", "scrutiny"),
                                   labels = c("Search", "Media", "Scrutiny"))

###!!!FIGURE 1.8!!!###
ggplot(data = df_mz_long_copy, aes(x = indicator, y = value, group = dyad)) +
  ggtitle("Mediatization") +
  geom_line(aes(linetype = Agreement, color = Agreement)) +
  geom_point(aes(shape = Agreement, color = Agreement)) +
  scale_shape_manual(values = c(1:28)) +
  scale_linetype_manual(values = c("solid", "dashed", "dotted", "dotdash", 
                                   "longdash", "twodash", "solid", "dashed", "dotted")) +
  scale_color_grey(start = 0.2, end = 0.8) +
  theme_minimal(base_family = "Cormorant", base_size = 14) +
  labs(x = "Indicator", y = "Value", color = "Agreement", linetype = "Agreement", shape = "Agreement") +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        legend.position = "right")

# Figure 1.9: Unclassified
###!!!FIGURE 1.9!!!###
create_trajectory_figure(unclassified, "Unclassified")

# Combine trajectory data
df_Ts <- rbind(unclassified, df_csd, df_es, df_pp, df_mz)

# =============================================================================
# CHAPTER 2: EXPLANATORY ANALYSIS (Figures 2.1-2.10)
# =============================================================================

# Load additional datasets
depth <- read.csv2(file="data/depth_details.csv") %>% rename(Agreement = agreement)
date <- read.csv2("data/date_data.csv") %>% rename(Agreement = agreement)
tci <- read.csv2(file="data/imports_relative.csv") %>% rename(Agreement = agreement)
values <- read.csv2(file="data/values_index.csv") %>% rename(Agreement = agreement)
search_results <- read.csv2(file="data/comm_power.csv") %>% rename(Agreement = agreement)

# Prepare explanatory datasets
deep_trade <- inner_join(df_Ts, depth)
parliamentary_assertion <- left_join(df_Ts, date)

# Globalization backlash data
tci <- subset(tci, Agreement != "brazil" & Agreement != "argentina" & 
              Agreement != "uruguay" & Agreement != "paraguay")
tci_long <- gather(tci, member_state, tci, at:uk, factor_key=TRUE)
tci_long <- na.omit(tci_long)
tci_long$dyad <- paste(tci_long$member_state, tci_long$Agreement)
tci_long$tci <- (tci_long$tci) * 100
base <- subset(df_Ts, select=c("dyad", "politicization", "Agreement", "trajectory"))
globalization_backlash <- left_join(tci_long, base)

# Existential threat data
values_long <- gather(values, member_state, inglehart_welzel, at:uk, factor_key=TRUE)
values_long$dyad <- paste(values_long$member_state, values_long$Agreement)
existential_threat <- left_join(values_long, base)

# Communicative power data
comm_power <- inner_join(df_Ts, search_results)

# =============================================================================
# FIGURES 2.1-2.2: DEEP TRADE HYPOTHESIS
# =============================================================================

# Figure 2.1: Level analysis
###!!!FIGURE 2.1!!!###
ggplot(deep_trade, aes(depth, politicization)) +
  geom_point(aes(shape = Agreement), size = 3, stroke = 1) + 
  geom_smooth(method = "lm", linetype = "dashed", color = "black", 
              fill = "grey70", alpha = 0.4) +
  scale_x_continuous(n.breaks = 10, limits = c(44, 92)) +
  scale_shape_manual(values = agreement_shapes, labels = custom_labels,
                     guide = guide_legend(ncol = 4, title = NULL)) +
  theme_minimal(base_family = "Cormorant", base_size = 14) +
  labs(x = "Regulatory Depth", y = "Level of Politicization") +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), legend.position = "top",
        legend.box.just = "center")

# Statistical analysis for Figure 2.1
model_depth_level <- lm(politicization ~ depth, data = deep_trade)
cat("Deep Trade Hypothesis - Level Analysis:\n")
cat("Estimate:", round(coef(model_depth_level)[2], 3), "\n")
cat("P-value:", round(summary(model_depth_level)$coefficients[2,4], 3), "\n")
cat("R-squared:", round(summary(model_depth_level)$r.squared, 3), "\n")
cat("95% CI:", round(confint(model_depth_level)[2,], 3), "\n\n")

# Figure 2.2: Trajectory analysis
deep_trade2 <- subset(deep_trade, trajectory != "unclassified")
###!!!FIGURE 2.2!!!###
ggplot(deep_trade2, aes(y = depth, x = reorder(trajectory, depth))) +
  geom_point(stat = "summary", fun = mean, size = 3) +
  geom_errorbar(stat = "summary", fun.data = mean_sdl, fun.args = list(mult = 1), width = 0.2) +
  theme_minimal(base_family = "Cormorant", base_size = 14) +
  labs(y = "Regulatory Depth", x = "Pattern") +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank())

# Statistical analysis for Figure 2.2
anova_depth <- aov(depth ~ trajectory, data = deep_trade2)
cat("Deep Trade Hypothesis - Trajectory Analysis:\n")
cat("F-statistic:", round(summary(anova_depth)[[1]][["F value"]][1], 3), "\n")
cat("P-value:", round(summary(anova_depth)[[1]][["Pr(>F)"]][1], 3), "\n")
if(summary(anova_depth)[[1]][["Pr(>F)"]][1] < 0.05) {
  cat("Tukey's HSD Test:\n")
  print(TukeyHSD(anova_depth))
}
cat("\n")

# =============================================================================
# FIGURES 2.3-2.4: PARLIAMENTARY ASSERTION HYPOTHESIS
# =============================================================================

# Figure 2.3: Level analysis
###!!!FIGURE 2.3!!!###
ggplot(parliamentary_assertion, aes(negotiation_start, politicization)) +
  geom_point(aes(shape = Agreement), size = 3, stroke = 1) + 
  geom_smooth(method = "lm", linetype = "dashed", color = "black", 
              fill = "grey70", alpha = 0.4) +
  scale_x_continuous(n.breaks = 13, limits = c(2007, 2020)) +
  scale_shape_manual(values = agreement_shapes, labels = custom_labels) +
  theme_minimal(base_family = "Cormorant", base_size = 14) +
  labs(x = "Start of negotiations", y = "Level of politicization", shape = NULL) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), legend.position = "top",
        legend.justification = "center", legend.box.just = "center") +
  guides(shape = guide_legend(ncol = 4))

# Statistical analysis for Figure 2.3
model_parl_level <- lm(politicization ~ negotiation_start, data = parliamentary_assertion)
cat("Parliamentary Assertion Hypothesis - Level Analysis:\n")
cat("Estimate:", round(coef(model_parl_level)[2], 3), "\n")
cat("P-value:", round(summary(model_parl_level)$coefficients[2,4], 3), "\n")
cat("R-squared:", round(summary(model_parl_level)$r.squared, 3), "\n")
cat("95% CI:", round(confint(model_parl_level)[2,], 3), "\n\n")

# Figure 2.4: Trajectory analysis
parliamentary_assertion2 <- subset(parliamentary_assertion, trajectory != "unclassified")
summary_data <- parliamentary_assertion2 %>%
  group_by(trajectory) %>%
  summarise(
    mean_negotiation_start = mean(negotiation_start, na.rm = TRUE),
    sd_negotiation_start = sd(negotiation_start, na.rm = TRUE)
  )

###!!!FIGURE 2.4!!!###
ggplot(summary_data, aes(x = reorder(trajectory, mean_negotiation_start), y = mean_negotiation_start)) +
  geom_point(size = 3) +
  geom_errorbar(aes(ymin = mean_negotiation_start - sd_negotiation_start, 
                    ymax = mean_negotiation_start + sd_negotiation_start), 
                width = 0.2) +
  theme_minimal(base_family = "Cormorant", base_size = 14) +
  labs(y = "Start of Negotiations", x = "Pattern") +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank())

# Statistical analysis for Figure 2.4
anova_parl <- aov(negotiation_start ~ trajectory, data = parliamentary_assertion2)
cat("Parliamentary Assertion Hypothesis - Trajectory Analysis:\n")
cat("F-statistic:", round(summary(anova_parl)[[1]][["F value"]][1], 3), "\n")
cat("P-value:", round(summary(anova_parl)[[1]][["Pr(>F)"]][1], 3), "\n")
if(summary(anova_parl)[[1]][["Pr(>F)"]][1] < 0.05) {
  cat("Tukey's HSD Test:\n")
  print(TukeyHSD(anova_parl))
}
cat("\n")

# =============================================================================
# FIGURES 2.5-2.6: GLOBALIZATION BACKLASH HYPOTHESIS
# =============================================================================

# Figure 2.5: Level analysis
###!!!FIGURE 2.5!!!###
ggplot(globalization_backlash, aes(tci, politicization)) +
  geom_point(aes(shape = Agreement), size = 3, stroke = 1) + 
  geom_smooth(method = "lm", linetype = "dashed", color = "black", 
              fill = "grey70", alpha = 0.4) +
  scale_shape_manual(values = agreement_shapes, labels = custom_labels,
                     guide = guide_legend(ncol = 4, title = NULL)) +
  theme_minimal(base_family = "Cormorant", base_size = 14) +
  labs(x = "Imports (relative to max value)", y = "Politicization") +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), legend.position = "top",
        legend.box.just = "center")

# Statistical analysis for Figure 2.5
model_glob_level <- lm(politicization ~ tci, data = globalization_backlash)
cat("Globalization Backlash Hypothesis - Level Analysis:\n")
cat("Estimate:", round(coef(model_glob_level)[2], 3), "\n")
cat("P-value:", round(summary(model_glob_level)$coefficients[2,4], 3), "\n")
cat("R-squared:", round(summary(model_glob_level)$r.squared, 3), "\n")
cat("95% CI:", round(confint(model_glob_level)[2,], 3), "\n\n")

# Figure 2.6: Trajectory analysis
globalization_backlash2 <- subset(globalization_backlash, trajectory != "unclassified")
###!!!FIGURE 2.6!!!###
ggplot(globalization_backlash2, aes(y = tci, x = reorder(trajectory, tci))) +
  geom_point(stat = "summary", fun = mean, size = 3) +
  geom_errorbar(stat = "summary", fun.data = mean_se, width = 0.2) +
  theme_minimal(base_family = "Cormorant", base_size = 14) +
  labs(y = "Imports (relative to max value)", x = "Pattern") +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank())

# Statistical analysis for Figure 2.6
anova_glob <- aov(tci ~ trajectory, data = globalization_backlash2)
cat("Globalization Backlash Hypothesis - Trajectory Analysis:\n")
cat("F-statistic:", round(summary(anova_glob)[[1]][["F value"]][1], 3), "\n")
cat("P-value:", round(summary(anova_glob)[[1]][["Pr(>F)"]][1], 3), "\n")
if(summary(anova_glob)[[1]][["Pr(>F)"]][1] < 0.05) {
  cat("Tukey's HSD Test:\n")
  print(TukeyHSD(anova_glob))
}
cat("\n")

# =============================================================================
# FIGURES 2.7-2.8: EXISTENTIAL THREAT HYPOTHESIS
# =============================================================================

# Figure 2.7: Level analysis
###!!!FIGURE 2.7!!!###
ggplot(existential_threat, aes(x = inglehart_welzel, y = politicization)) +
  geom_point(aes(shape = Agreement), size = 3, stroke = 1) + 
  geom_smooth(method = "lm", linetype = "dashed", color = "black", 
              fill = "grey70", alpha = 0.4) +
  scale_shape_manual(values = agreement_shapes, labels = custom_labels,
                     guide = guide_legend(ncol = 4, title = NULL)) +
  theme_minimal(base_family = "Cormorant", base_size = 14) +
  labs(x = "Difference on survival/self-expression index", y = "Level of politicization") +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), legend.position = "top",
        legend.box.just = "center")

# Statistical analysis for Figure 2.7
model_exist_level <- lm(politicization ~ inglehart_welzel, data = existential_threat)
cat("Existential Threat Hypothesis - Level Analysis:\n")
cat("Estimate:", round(coef(model_exist_level)[2], 3), "\n")
cat("P-value:", round(summary(model_exist_level)$coefficients[2,4], 3), "\n")
cat("R-squared:", round(summary(model_exist_level)$r.squared, 3), "\n")
cat("95% CI:", round(confint(model_exist_level)[2,], 3), "\n\n")

# Figure 2.8: Trajectory analysis
existential_threat2 <- subset(existential_threat, trajectory != "unclassified")
###!!!FIGURE 2.8!!!###
ggplot(existential_threat2, aes(x = reorder(trajectory, inglehart_welzel), y = inglehart_welzel)) +
  geom_point(stat = "summary", fun = mean, size = 3, color = "black") +
  stat_summary(fun.data = mean_cl_normal, geom = "errorbar", width = 0.2, color = "black") +
  theme_minimal(base_family = "Cormorant", base_size = 14) +
  labs(y = "Difference on survival/self-expression index", x = "Pattern") +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank())

# Statistical analysis for Figure 2.8
anova_exist <- aov(inglehart_welzel ~ trajectory, data = existential_threat2)
cat("Existential Threat Hypothesis - Trajectory Analysis:\n")
cat("F-statistic:", round(summary(anova_exist)[[1]][["F value"]][1], 3), "\n")
cat("P-value:", round(summary(anova_exist)[[1]][["Pr(>F)"]][1], 3), "\n")
if(summary(anova_exist)[[1]][["Pr(>F)"]][1] < 0.05) {
  cat("Tukey's HSD Test:\n")
  print(TukeyHSD(anova_exist))
}
cat("\n")

# =============================================================================
# FIGURES 2.9-2.10: COMMUNICATIVE POWER HYPOTHESIS
# =============================================================================

# Figure 2.9: Level analysis
###!!!FIGURE 2.9!!!###
ggplot(comm_power, aes(bilaterals, politicization)) +
  geom_point(aes(shape = Agreement), size = 3, stroke = 1) + 
  geom_smooth(method = "lm", linetype = "dashed", color = "black") +
  scale_x_continuous(n.breaks = 10, limits = c(0, 1900)) +
  scale_shape_manual(values = agreement_shapes, labels = custom_labels,
                     guide = guide_legend(ncol = 4)) +
  theme_minimal(base_family = "Cormorant", base_size = 14) +
  labs(x = "Returns on Bilaterals.org", y = "Level of politicization") +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), legend.position = "top",
        legend.title = element_blank(), legend.box.just = "center")

# Statistical analysis for Figure 2.9
model_comm_level <- lm(politicization ~ bilaterals, data = comm_power)
cat("Communicative Power Hypothesis - Level Analysis:\n")
cat("Estimate:", round(coef(model_comm_level)[2], 3), "\n")
cat("P-value:", round(summary(model_comm_level)$coefficients[2,4], 3), "\n")
cat("R-squared:", round(summary(model_comm_level)$r.squared, 3), "\n")
cat("95% CI:", round(confint(model_comm_level)[2,], 3), "\n\n")

# Figure 2.10: Trajectory analysis
comm_power_t <- subset(comm_power, trajectory != "unclassified")
###!!!FIGURE 2.10!!!###
ggplot(comm_power_t, aes(x = reorder(trajectory, bilaterals), y = bilaterals)) +
  stat_summary(fun = mean, geom = "point", shape = 20, size = 3, color = "black") +
  stat_summary(fun.data = mean_cl_normal, geom = "errorbar", width = 0.2, color = "black") +
  theme_minimal(base_family = "Cormorant", base_size = 14) +
  labs(x = "Pattern", y = "Returns on Bilaterals.org") +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank())

# Statistical analysis for Figure 2.10
anova_comm <- aov(bilaterals ~ trajectory, data = comm_power_t)
cat("Communicative Power Hypothesis - Trajectory Analysis:\n")
cat("F-statistic:", round(summary(anova_comm)[[1]][["F value"]][1], 3), "\n")
cat("P-value:", round(summary(anova_comm)[[1]][["Pr(>F)"]][1], 3), "\n")
if(summary(anova_comm)[[1]][["Pr(>F)"]][1] < 0.05) {
  cat("Tukey's HSD Test:\n")
  print(TukeyHSD(anova_comm))
}
cat("\n")

# =============================================================================
# MULTIPLE REGRESSION ANALYSIS
# =============================================================================

# Prepare combined dataset
test <- df_Ts %>% dplyr::select(dyad, politicization, trajectory)
parliamentary_assertion_test <- parliamentary_assertion %>% dplyr::select(dyad, negotiation_start)
deep_trade_test <- deep_trade %>% dplyr::select(dyad, depth)
globalization_backlash_test <- globalization_backlash %>% dplyr::select(dyad, tci)
existential_threat_test <- existential_threat %>% dplyr::select(dyad, inglehart_welzel)
comm_power_test <- comm_power %>% dplyr::select(dyad, bilaterals)

test <- inner_join(test, parliamentary_assertion_test, by="dyad") %>%
  inner_join(deep_trade_test) %>%
  inner_join(globalization_backlash_test) %>%
  inner_join(existential_threat_test) %>%
  inner_join(comm_power_test)

# Rename variables
test <- test %>% rename(
  imports = tci, framing = bilaterals, values = inglehart_welzel,
  inst_change = negotiation_start
)

# Mean center predictors
mean_center <- function(df) {
  df_centered <- as.data.frame(scale(df, center = TRUE, scale = FALSE))
  return(df_centered)
}

predictors <- test[, c("inst_change", "depth", "imports", "values", "framing")]
outcome <- test$politicization
predictors_centered <- mean_center(predictors)
data_centered <- cbind(predictors_centered, outcome)

# Table 2.1: Multiple Regression Models
model_1 <- lm(outcome ~ inst_change + depth + imports + values + framing, data = data_centered)
model_2 <- lm(outcome ~ depth + imports + values + framing, data = data_centered)
model_3 <- lm(outcome ~ depth * framing + imports * framing + values * framing + framing, data = data_centered)

###!!!TABLE 2.1!!!###
stargazer(model_1, model_2, model_3, type = "text", title = "Multiple Regression Results")

# VIF Analysis
vif_full <- vif(model_1)
vif_reduced <- vif(model_2)

cat("Variance Inflation Factors:\n")
cat("Full Model:\n")
print(round(vif_full, 3))
cat("\nReduced Model:\n")
print(round(vif_reduced, 3))

# =============================================================================
# MULTINOMIAL LOGISTIC REGRESSION (Ridge Regression)
# =============================================================================

# Prepare data for multinomial analysis
test_filtered <- subset(test, trajectory != "unclassified")
test_filtered$depth_centered <- scale(test_filtered$depth, center = TRUE, scale = FALSE)
test_filtered$imports_centered <- scale(test_filtered$imports, center = TRUE, scale = FALSE)
test_filtered$values_centered <- scale(test_filtered$values, center = TRUE, scale = FALSE)
test_filtered$framing_centered <- scale(test_filtered$framing, center = TRUE, scale = FALSE)

# Residualize framing with respect to values to address multicollinearity
framing_resid <- residuals(lm(framing_centered ~ values_centered, data = test_filtered))
test_filtered$framing_resid <- framing_resid
test_filtered$pattern <- as.factor(test_filtered$trajectory)

# Prepare data for ridge regression
x <- as.matrix(test_filtered[, c("depth_centered", "framing_resid", "imports_centered", "values_centered")])
y <- as.numeric(test_filtered$pattern)

# Run ridge regression with cross-validation
class_weights <- table(test_filtered$pattern)
ridge_model <- cv.glmnet(x, y, alpha = 0, family = "multinomial", weights = 1/class_weights[y])

# Extract coefficients
coef_ridge <- coef(ridge_model, s = "lambda.min")
print(coef_ridge)

# Predict using the ridge model
predictions <- predict(ridge_model, newx = x, s = "lambda.min", type = "class")

# Create confusion matrix
confusion_matrix <- table(Predicted = predictions, Actual = test_filtered$pattern)

###!!!TABLE 2.4!!!###
cat("Ridge Regression - Confusion Matrix:\n")
cat("=====================================\n")
print(confusion_matrix)

cat("\nAnalysis completed. All essential figures (1.1-1.9, 2.1-2.10) and tables (2.1, 2.4) have been generated.")